Markov model for multispectral image analysis: application to small magellanic cloud segmentation

This paper deals with the unsupervised segmentation of astronomical multiband images. Most of these images have the particularity to be quantized on float numbers with large luminance range, on different wavelengths. These characteristics require to manipulate large amount of extremely accurate data on each spectral band, which is very different to the case of 8-bits-integer coded pixels. We present some results obtained on Multiwavelength images of the Small Magellanic Cloud, by using the Marginal Posterior Mode (MPM) estimator on a quadtree structure under Markovian assumption. The estimation of the model parameters is then addressed with Expectation-Maximization (EM)-type algorithms, allowing unsupervised hyperparameter estimation. The main interest of this modeling effort lies in its generality : the algorithm handles multiwavelength floating data in a single upward and downward scan on the quadtree. A new aspect in this paper concerns the noise statistics that are supposed to be lognormal for each class. Another new aspect, is the in-scale-coding of the label map. 1. MULTIWAVELENGTH IMAGES IN ASTRONOMY The study of star formation mechanisms and their relationship with interstellar medium is one of the most dynamic fields of research in astronomy. The interstellar medium is composed of a mixture of gas (mainly hydrogen and helium) in different phases and dust grains (mainly carbon and silicium). Stars are forming from the gas in the heart of molecular Hydrogen complexes. The molecular clouds themselves result from the formation of molecules via chemical reactions in dense atomic gas. The ultra violet emission of massive newborn stars ionizes the surrounding gas, giving birth to the so-called HII regions. Observations at various wavelengths are necessary to study the different states of gas and dust and their mutual relationships as well as their links to star formation. The emission line at 21cm in the radio wavelength range gives the column density of atomic hydrogen (HI), which is the number of hydrogen atoms in a unitary section cylinder along the line of sight. Emissions in the far Infrared at 100 or 170 microns are directly due to the thermic emission of big dust grains. The Hα emission, at 656.3nm, is the strongest Hydrogen recombination line in the HII regions and is generally superimposed on the optical continuum of stars[1]. Molecular Hydrogen are difficult to detect, but estimation of their density can be infered from the intensity of the CO lines at millimeter wavelengths, assuming that this CO and molecular hydrogen are formed at quite the same time. The images studied here are centered at equatorial coordinates : 00h52m07s -72d49m (J2000) covering a large fraction of the SMC. These are the following : • The ”HI image”, Figure 3a, observed at 21.1cm is obtained from both the Australia Telescope Compact Array (ATCA) radio interferometer and from the Parkes telescope. These combined observations allow a good spatial frequency coverage ; • The ”PHOT image”, Figure 3b, coding the luminosity value recorded at 170μm by the ISOPHOT photometer on board ISO, the Infrared Space Observatory operated by the European Space Agency (ESA). On this particular image, some missing pixels appear as a regular lattice of undefined pixels values. These gaps come from accidental undersampling of the observations ; • The ”Hα image”, Figure 3c, observed at 656.3nm for the same region is used as explanatory data in order to examine correlation with the star forming regions ; • The ”IRAS100 image”, observed at 100μm from the IRAS instrument (Infrared Astronomical Satellite) is used in conjunction with the PHOT image to estimate the temperature of regions of different components. All these images have been registered and resampled with respect to the ISOPHOT image. The goal consists in using ”HI image” and ”PHOT image” to segment the SMC cloud into different classes, in order to be able to estimate physical parameters such as the temperature for each of them[2]. In order to allow that, we use a bi-wavelength Markovian-in scale segmentation in order to generate a segmentation map labelling the pixel in the different classes. 2. IMAGE ANALYSIS WITH THE MARKOVIAN HIERARCHICAL METHOD For segmentation tasks, with strongly noisy images, Markovian assumption allows the description of global behaviors by considering, on a predetermined spatial neighboorhoud around each pixel, statistical relationships between observation field and label field (i.e., segmented image). Nevertheless, Markov Random Fields (MRFs) lead to robust but iterative procedures (i.e., computionally demanding due to slow convergence[3]). To avoid such difficulties, hierarchical modeling allows the definition of different coarse-to-fine strategies, under (spatially or/and in-scale) Markovian assumptions[4]. In this way, we consider a special class of Markov models, which helps to circumvent the latter drawbacks (i.e., iterative often computational intensive estimation algorithms). Indeed, the Causal-in-scale Markov Random models, attached to the nodes of a quadtree, are defined in a general manner, allowing to analyse simultaneously multiband and/or multiresolution images [5]. The great advantage of this model results in an interesting causality property through scale, which allows the design of exact and non-iterative inference algorithms which are similar to those used in the context of Markov Chain Models (MCM). Indeed, another possible approach consists in scanning the image with a fractal path (scan of Hilbert-Peano for example) : in the case of multiband pictures, the chain is composed of a succession of multiwavelength observations on which a Markovian chain model is applied[6]. Nevertheless, the Markov Chain approach is not suited in the case of missing observations on the sampled grid (Figure 2), whereas hierarchical approach on the quadtree can deal with such unobserved samples[2]. Indeed, the Markovian quadtree segmentation has the capability to overcome missing zones on the observations (data likelihood equals one for such pixels) and even to provide a segmentation map without missing labels : the quadtree propagation of the labels and transitions laws in the method, allows to overcome the missing pixels problem and spacing artefacts, as shown on ”PHOT image”, Figure 3(b). In our example, on both HI and ISOPHOT images, we applied a markovian quadtree segmentation. A new aspect in this paper concerns the noise statistics that are supposed to be lognormal for each class. The data-driven model follows a lognormal probability density function, whose parameters are estimated by using an ICE procedure[7] using the Maximum Likelihood estimators which are presented in the next section. We obtain the 8-classes segmentation map given in Figure 4(a). Another new aspect, is the coding of the label map. In fact, we code the obtained segmentation map at s = 8 different scales for a maximum of 16 classes (4 bits) for each, from the finest resolution s = 0. Thus, each site receives a label coded on 32 bits : 4 bits per scale on 8 scale levels. The observation of the different map in scale is of great interest in the context of astronomy where a single site should belong to a class which depends of the scale level : for example a star at full resolution belongs to a galaxy at a coarser resolution and so on... multi-scale label on each site of the grid (32 bits) z }| { class 5 coarser s=7 z}|{ 0101 − class 7 s=6 z}|{ 0111 − class 6 s=5 z}|{ 0110 −...− finest s=0 class 4 z}|{ 0100 Figure 5(a-d) present the label maps for four different coarser scales. All the regions identified here are homogeneous with respect to their spectrophotometrical information. Futher analysis of physical properties in the classes can then take place, particularly the class by class correlated temperature assumption with the IRAS map[2]. 3. LOGNORMAL PDF FOR THE DATA-DRIVEN TERM The lognormal distribution fY is used extensively in reliability applications to model failure times and in astronomy to model noisy observations y for a given class : the formula below are with location parameter equals to a, scale parameter equals to m and shape parameter equals to σ : fY (y) = 1 (y − a) σ√2π exp −(ln(y−a)−m)2 2σ2 ; y > a with E [y] = exp σ2 2 +a and V ar [y] = exp 2 (exp 2 −1) A variable Y = exp + a is lognormally distributed if X is normally distributed with ’ln’ denoting the natural logarithm. The maximum likelihood estimates for the scale parameter m, the shape parameter σ and the location parameter a, are b m = PN i=1 ln (yi) N and b σ = vutPNi=1[ln(yi)− PN j=1 ln yj N ]2 N − 1 ba = min (yi)∀i∈[1,N ] If the location parameter is known, it can be subtracted from the original data points before computing the maximum likelihood estimates of the shape and scale parameters. The segmentation maps Figure 5 are generated with such a data-driven model, which fits with a large variety of distributions as shown in Figure 1. 4. CONCLUSION AND PERSPECTIVES We can see on Figure 5 the label map obtained for different scales : this allows an in-scale labelling process and open the way to the definition of a vector of labels for each site on the grid, corresponding to the class associated to each multiwavelength pixels, at different scales of resolution. This approach is of major interest in astronomy, where each pixel has to be be labelled according to the scale of observation. Thus, by combining the segmentation map and the missing pixels list, we derive a tool for selecting the valid pixels of the different classes. For each of them, statistical, geometrical or physical parameters can be evaluated, at different scale. Acknowledgments: The authors thank C. Bot and F. Bonnarel for fruitfull discussions, (CDS Strasbourg Astronomical Observatory) about astronomical image interpretation and management. This research is supported by the French government (ACI-Grid / IDHA project: Action Concertée Incitative-Globalisation des Ressources Informatiques et des Données / Images Distribuées Hétérogènes pour l’Astronomie, 2001-2003) 5. REFERENCES [1] C. Bot, F. Boulanger, K. Okumura, and B. Stepnik, “Multiwavelength analysis of the dust emission in the small magellanic cloud,” Proc. of symp. ”Exploiting the ISO data archive Internet Astronomy in the Internet Age, vol. C. Gry et al. eds., ESA SP 511, 24-27 June 2002. [2] M. Louys, A. Oberto, C. Bot, and C. Collet, “Hierarchical markovian inference for astronomical multiband images segmentation,” Physics in Signal and Image Processing, PSIP 2003, vol. 1, pp. 61–64, January 2003, Grenoble, France. [3] S. Geman and D. Geman, “Stochastic relaxation, Gibbs distributions and the Bayesian restoration of images,” IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. PAMI-6, no. 6, pp. 721–741, November 1984. [4] P. Pérez and F. Heitz, “Restriction of a markov random field on a graph and multiresolution statistical image modeling,” IEEE Trans. Information Theory, vol. 42:1, pp. 180–190, 1996. [5] J.-M. Laferté, P. Pérez, and F. Heitz, “Discrete markov image modeling and inference on the quad-tree,” IEEE Trans. Image Process., vol. 9, no. 3, pp. 390–404, March 2000. [6] N. Giordana and W. Pieczynski, “Estimation of generalized multisensor hidden Markov chains and unsupervised image segmentation,” IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 19, no. 5, pp. 465–475, 1997. [7] W. Pieczynski, “Statistical image segmentation,” Machine Graphics and Vision, vol. 1, no. 2, pp. 261–268, 1992. 0 2 4 6 8 10 12 14 16 18 0 0.2 0.4 0.6 0.8 1 1.2 1.4 Lognormal distribution